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Abstract 

In canonical quantum cosmology, the wave function of the universe lacks explicit time depen¬ 
dence. However, time evolution may be present implicitly through the semiclassical superspace 
variables, which themselves depend on time in classical dynamics. In this paper, we apply this 
approach to an oscillating universe model recently introduced by Graham et al. By extending the 
model to include a massless, minimally coupled scalar held f which has little effect on the dynamics 
but can play the role of a “clock”, we determine the decay rate of the oscillating universe. 
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I. INTRODUCTION 


Oscillating cosmological models have been extensively studied over the years (see, e.g., 
and references therein). Such models can be classically stable with respect to small 
perturbations M, but it was pointed out in [5], [6] that they are generically unstable with 
respect to decay through quantum tunneling to zero size. Our goal in the present paper is 
to calculate the corresponding decay rate. 

In the semiclassical regime, the decay rate can be expressed as 

r = ^ e - 2|5£l , (i) 

where Se is the under-barrier Euclidean action. The action Se has been calculated in 
Refs. J5, E] for some simple FRW models. Here, we would like to go beyond that and also 
calculate the pre-exponential factor A, at least in the framework of the FRW models under 
consideration. 

The problem we have to address is that the rate T is the decay probability per unit time, 
and the time variable is conspicuously absent in the formalism of quantum cosmology. In 
particular, the wave function of the universe T is independent of time. The way out of 
this impasse was already pointed out by DeWitt [7], who suggested that the role of a clock 
could be played by some geometric or matter variable. Any time evolution should then 
be understood implicitly, in terms of the canonical variables themselves. We shall adopt 
this approach here and use it to calculate the decay rate in the “simple harmonic universe” 
(SHU) model of Ref. [2], suitably extended to include a clock. 

The SHU is a closed universe with energy density due to a negative cosmological constant 
A < 0 and a matter component having equation of state P = wp. Oscillating solutions exist 
for — 1 < w < —1/3. The gravity of matter with such equation of state is repulsive; it causes 
a contracting universe to bounce and to start expanding. On the other hand, a negative A 
causes an expanding universe to turn around and start contracting. As a result, the universe 
will oscillate between some maximum and minimum size. Here, we focus on the case where 
w = —2/3 for calculation simplicity, though the situation is qualitatively similar for other 
values of w. 

The FRW model of Ref. [2] includes a single dynamical variable - the scale factor a. 
In the case of an oscillating universe, the scale factor evolution is not monotonic; hence it 
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cannot serve as a time variable. We therefore introduce a second minisuperspace variable 
- a homogeneous, massless, minimally coupled scalar field 0, which will play the role of a 
clock. We shall assume that the contribution of this field to the total energy density of the 
SHU is negligible, so that its presence has little effect on the dynamics of oscillations and 
does not alter the stability analysis of [3]. 

The rest of the paper is organized as follows. In the next section we review the SHU model 
and extend it by introducing a homogeneous scalar field. The semiclassical wave function 
for this model is found in Sections III and IV. The tunneling rate and the corresponding 
lifetime of the universe are calculated in Sec. V. This rate has the expected relation to the 
tunneling probability, as we show in Sec. VI. Finally, in Sec. VII we give a brief summary 
and some concluding remarks. 


II. SIMPLE HARMONIC UNIVERSE WITH A CLOCK 


A. Simple harmonic universe 


Classical dynamics of the SHU are described by the k = +1 Friedmann equation 

a 2 1 8itG 


a 2 a 2 


-P(a), 


with energy density 


a 


p(a ) = A + -. 

a 


( 2 ) 


( 3 ) 


Here, A < 0 is the cosmological constant and the second term in (|3j) describes “domain wall 
matter” with equation of state w = —2/3. We shall assume that both parameters A and a 
are small in Planck units, 

|A| < G~ 2 , a < G~ 3/2 . (4) 


The Friedmann equation has oscillating solutions for the scale factor 


i(t) — oj x (y — y/q 2 — 1 cos(c ot)), 


( 5 ) 


where 


to = 


87 tG 


|A| 


7 = 


' 2i tG a 2 

TjAf 


( 6 ) 


3 








and we have to require that 7 > 1. The universe oscillates between maximum and minimum 
values of 


a±=u 1 (7 ± vV - l) • (7) 

Alternatively, the system also may be described by the constrained Hamiltonian 

(8) 

where the momentum conjugate to a is 

371 ■ ( q^ 

Pa = ~2 G aa 

and the potential U (a) is 

U(a) = a 2 ^1 - ^a 2 p(a)^j . ( 10 ) 

In the context of quantum cosmology, the momentum becomes the differential operator 
p a —y —id/da and the wave function of the universe T(a) satisfies the Wheeler-deWitt 
equation 

HV(a) = 0 . ( 11 ) 


B. SHU with a scalar field 


We now modify the SHU model by adding a homogeneous, massless, minimally coupled 
scalar field pit). The corresponding Hamiltonian constraint is 

G ( „ 3 


n = 


3tt a 


pI 


4ttGg 2 ^ 


Pu> + U(a) \ — 0 , 


where 


P<p = 2tt 2 o?p 


( 12 ) 


(13) 


is the momentum conjugate to 0. The momentum p$ is a constant of motion, p^ — 0. 


Without loss of generality, we shall assume that p v > 0. Then it follows from Eq. (13) that 
the scalar field ip increases monotonically; hence it can be used as a time variable. 

In quantum cosmology, we make the replacement p a —> —id/da and p v —> —id/dcp, and 
the Hamiltonian constraint H = 0 becomes the Wheeler DeWitt (WDW) equation 

3 d 2 


+ “ 2f/(a) + 4 xGV, 


T(a, 0) = 0. 


(14) 
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Here, we have adopted the ordering of the non-commuting factors a and d/da proposed in 


[ 8 ], for which the differential operator in Eq. (14) becomes a covariant Laplacian. 

With the change of variables a = In (orya), 0 = ( 47 tG/ 3 ) 1//2 <^, the WDW equation becomes 


92 TTf ^ , 92 

Ufa) + 


da 2 


= 0, 


where the potential a 2 U(a ) = 17(a) is 

U{a) = /T 2 e 4a (1 - 2 e Q + 7 “ 2 e 2 “) , 

where 

/ 2G\ 2 2 327T 3 2 

B = — uW = —-G A o l . 

A V 3vr y 1 27 


(15) 


(16) 


(17) 


The WDW equation (15) separates, and the general solution can be expresses as a su¬ 


perposition of terms of the form 


*= e' r *f r (a), 


(18) 


where the separation parameter p is the eigenvalue of the momentum and the function 
f p (a) satisfies the equation 


d 2 


+ U p {a) 


fp(<*) = 0 


with 


da 2 

U p (a) = U(a) — p 2 . 


(19) 

( 20 ) 


The effective potential U p (a) is plotted in Fig. 1. 

We see that inclusion of a scalar held has the effect of decreasing the potential U (a) by a 
constant term, — p 2 . We assume that this term is small compared to the characteristic scale 
of the potential, that is, 

p^fi- 1 . ( 21 ) 

This term, however, does have an effect near the turning points 

a± = In (q 2 ± 7 2 y/l - q" 2 j ( 22 ) 

of the unperturbed potential, U(a±) = 0. The turning points in the presence of a scalar 
held, 


«i = a_ — 5a\ 
c>f 2 = a + T Sa 2 
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(23) 

(24) 













can be found by solving U p (a) = 0. To the lowest order in p 2 , we have 


Sa i — — 

1 l TV 


P 


\U'(a-)\ 


5a 2 ~ — 


P 


U'(a+Y 


U\a _) = —— (-1 + y/l - 7 2 ) 7 8 (1 + (-1 + \/l - 7 2 ) 7 2 


The derivatives of the potential appearing in Eqs. (|25|),(26) are 

1 / ,-x4 

J 2 

U\a + ) = + y/l -7" 2 ) 7 8 (-1 + (i + V 1 -7 -2 ) 7 2 ) , 

or, by order of magnitude, 


U\a_) ~ -^“ 2 , C/'(a+) ~ f3~ 2 7 10 . 


(25) 

(26) 

(27) 

(28) 

(29) 


Since 7 > 1, we can write 

5a± < /3 2 p 2 -C 1. (30) 

Apart from shifting the turning points a±, the scalar held also modifies the character of 
the potential at small a (a —> — 00 ), introducing another classically allowed region (region 1 
in Fig. 1). The potential at a —> —00 can be approximated as U p (a) ~ f3~ 2 e Aa — p 2 , so the 
boundary of this region is approximately 

«o-^ln {Pp). (31) 

In order to justify semi classical treatment, we shall require that the corresponding scale 
factor is large in Planck units, 

00 = (^) 1/2 »C? 1/2 . (32) 

which implies p^> 1 . 

The two classically allowed regions are separated by a barrier (region II) extending be¬ 
tween the turning points cxo and aq. The situation is therefore analogous to a particle in a 
metastable state: the particle is localized in an approximate energy eigenstate, but the cor¬ 
responding energy eigenvalue takes on an imaginary part, indicating a non-vanishing decay 
rate. We shall see that the wave function in our model exhibits a very similar behavior. 
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y = 1.03, (3 = .03, p = 3 


U p (a) 



FIG. 1. Effective potential U p (a), with 7 = 1.03, f3 = .03, and p = 3. 


III. SEMICLASSICAL SOLUTIONS 


Solutions to the WDW equation are specified by imposing boundary conditions appro¬ 
priate to the problem. Here, we first require that the wave function must decay under the 
infinite barrier to the right of 0 : 2 , / p («) —» 0 as a —* 00 (this is the same boundary condition 
chosen in 0 ). In addition, we require that the solution be outgoing (left-moving) in the 
region a < a 0 . Physically, this means that once the oscillating universe tunnels through the 
barrier, it collapses to the singularity at a = 0 (a —>• — 00 ). In other words, the singularity 
is a point of no return: the probability of getting back from a = 0 to a = cio is zero. 

Sufficiently far from the turning points, we can determine the solutions using the semi- 
classical approximation: 


/„(<,) ~ . Cie , "(b e+y.. + C2e ‘* /1 


\-U r (a)Y ' 1 


[-C,(a )] 1/4 


yJ-U p {a')da' 


(33) 


for left (+) and right (—) moving waves in the classically allowed region, U p (a ) < 0, and 


/p(«) - 


D 1 


(W ) 1/4 


V U ^ da ' + 


Do 


(U p (a))V* 




for growing (+) and decaying (—) solutions under the barrier, U p (a) > 0. 


( 34 ) 
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Near a turning point, a = a*, wher^] 


Re U p (a *) = 0, 


(35) 


the semiclassical approximation breaks down. In such regions we use the standard technique 
[9] and approximate the potential by a linear function, 


U(a * + 5a) ~ U(a*) + U'(a*)5a = U'(a*)(a — a*). 


(36) 


Setting z = ([/'(a;*)) 1 / 3 (ck — a*), the approximate WDW equation near a turning point is 



(37) 


The solution is a linear combination of Airy functions Ai(z) and Bi(z), having the asymptotic 
{z —* oo) forms 



(38) 


(39) 


Ai(-z) ~ ^-z-^shA-z*/ 2 + - 

V 71 L 3 4 . 

Bi(-z) ~cos [-^ 3/2 + - 
A /7T 3 4 


(40) 


(41) 


On the other hand, for the linearized potential, 


■,) 1/2 da= 2 l t/ '(° - )] 1/2 ( a - a.) 3 / 2 
o 


[Cfa)] 1 ^ 2 ^ a; [JUa,)] 1 / 2 (a 


(42) 



= A" 


(43) 


We thus see that the linearized approximation near the turning points, with Airy function 
solutions, matches onto the WKB solutions away from the turning points. We determine 
solutions in all regions by imposing the boundary conditions in regions I and IV, and match 
semiclassical solutions across the turning points ao, a± and a^-- 

We first apply the boundary condition in region IV, to the right of CK 2 . There, the solution 
consists only of a decaying mode, f p (a —» oo) —» 0: 



(44) 


1 We have to use the real part, since the parameter p is generally complex, and thus the potential U p (a ) is 
also complex. We shall see, however, that the imaginary part of p is exponentially suppressed. Hence we 
shall disregard it everywhere, except for the calculation of the decay rate. 











where A = const. With the asymptotic form of the Airy functions, this fixes the coefficients 
across o: 2 in region III: 



e ~ i7r / 4 A 

2[-t/ p (a)]i/4 


y/-Up(a')da' _j_ ^-i f “ 2 y/-U p (a')da' j 


(45) 


The second boundary condition - that the solution be outgoing in the region a < a 0 
means that the solution must take the form 


/p(«) = 


e ^ y/—Up(a l )da l 


(46) 


-£/»] 1/4 ^ 

where I? is a constant coefficient 0 We can now use the same method as above to match the 
solutions across the turning point «o and fix the coefficients of the wave function in region 
II: 


/"(«) = 


-iB 


2 [ ty «)] i/4 
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[[y p (a)] 1/4 


X 0 \/wW 


(47) 


We now have expressions for solutions f J p and f^ 1 in terms of coefficient B, and solutions 
f p u and f p l in terms of coefficient A ; we must now reconcile solutions everywhere in terms 
of a single coefficient. The general solution to the left of op is 


/"(«) = 


A' 


,-fa 1 VUpWM + 


5' 


J7 1 A /u v (a')da' 


(48) 


2 [ U P { a ) Y /A [^(«)] 1/4 

Matching this across aq with the aid of the linearized approximation, we find the form of 
the solution to the right of an, 

3 —27t/4 


fp » = 


((A' + wyK i V- u A*')d«' + + / 


—U p (a')da' 


(49) 


2[-U r (a)Y!‘ 

The coefficients A 1 and B' can now be determined by noting that the solution f p r (a) 
in Eq. (48) must match the solution in Eq. (47) determined with the boundary conditions. 

K = I \/u p (a')da', (50) 


Defining 


' «0 


we find the relations 


A' = 2 Be K 


B' = 



(51) 

(52) 


2 Note that two linearly independent solutions in the limit a —> —oo are f p (a) oc exp(±ipa), so the outgoing 
mode can be unambiguously identified. 
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Then the solution in region III is 


fp («) = 




2[~U P (a)} 1/4 


(^2e K + e <; “i V~ u P W da ' + i ^2e K - e~ if ^ V~ u »W da ' 

(53) 


Similarly, we require that the solutions fl H (o>) from Eq. (45) and Eq. (53) agree: 


B 


((2e* + <q£) e ‘A + i ( 2 e* - e -<C, (54) 

= A y/- u p( a ') dOL ' + V z uAA) da '^J . (55) 


Defining 


the relations 


J 



\J —U p (a')da', 


A 

A 


iB 




-iB 



(56) 


(57) 

(58) 


must be simultaneously satisfied. 

In order for a solution to exist, we must require 


or 


1 + 


= — e 


2 ij 



(59) 


(60) 


where n is an integer. When K is large, expanding the logarithm results in the approximate 
relation 


J ~ 



7r + -e 
4 


-2 K 


(61) 


This condition will later be used to determine the momentum eigenvalue p and the decay 
rate T. 
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IV. EVALUATION OF J AND K 


In the semiclassical regime, assuming that 7 is not too close to 1 and that p ^ (3 1 , the 
contribution to J from p may be treated as a perturbation: 


J — / \J—U (cc) + p 2 da 

' a± 

r&+ _ 

\/—U (a) da + 


(62) 


p 


, - —dcx ~f~ 

_ 2 \J—U (a) \j ai 


ca 2 


= Jq + J\ + J‘2 ■ 

We may then evaluate J 0 , Ji, and J 2 analytically: 

Jo = 


_ p* 
Jl ~ 2 


\/—U (a) da + j \J— U (a)dajQ?>) 

(64) 


J‘2 ~ 


\J—U(a)da = ^(7 3 - 67 s + 57 ? ), 

( 65 ) 

r 1 i w 

/ —. da = -, 

L. v / - f/ («) 2 

(66) 

pot 2 

y/—U(a)da+ / y/—U(a)da 

( 67 ) 

J a+ 


V-f/'(a_) 5 a? /2 + ^/U'(a + )Sa 2 /2 ^ 

(68) 

P3 + p3 ^ ~P 2 p 3 . 

J( 7 '(a_)| f 7 '(a + )J PP 

( 69 ) 


Here, in the calculation of J 2 we have expanded t/(a) near ck± and used Eqs. (25), (26) and 


p9j). The contribution from the correction to the turning points, J 2 , is small compared to 

7 Tf3p 2 


Ji, hence we can write 


J - Jo + 


(70) 


To evaluate K, we again expand about the p = 0 limit, K = K 0 + 8K, where 

1 


K n = 


\/U(a)da = 


24 (3 


157® - 137* + ? ( 7 3 - 67 s + 57 7 ) hi (^-T 


(71) 


and 5K ~ (3p 2 includes all corrections due to p. Even though 5K <C Ji 0 , we cannot generally 
neglect 8K. Since p 3> 1 is required for our semiclassical analysis, we can have 5K > 1 even 


if (3p -C 1. Neglecting SK in Eq. (61) is justified only if 

I3p 2 « 1, 


( 72 ) 


To simplify further analysis, we shall assume this condition to be satisfied. 
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V. THE DECAY RATE 


We shall now use Eqs. (61) and (70) to determine the momentum eigenvalue p. We first 
assign to p real and imaginary parts, 


p = p + ip , 


(73) 


with p' and p" real. We shall assume that p" <C p'\ this will be justified below. (Note that 
we also neglected the effect of p" on the classical turning points in earlier sections.) 


Substituting J from (70) in (61), using (73) and neglecting p " 2 compared top' 2 , we obtain 
two relations: 


W 2 f 1 

^ 2 = f n ^ 2 _ X 

TT0p'p" = le- 2K ». 


With J 0 from Eq. (65), the first of these relations becomes 


7 3 (q 2 — l)(5y 2 — 1) = 4f3[(2n + 1) — f3p n ] a 4/3(2n + 1). 


(74) 

(75) 


(76) 


Disregarding the small correction introduced by the “clock”, as we did in the last step, 
this is a quantization condition on the parameters of the model /3 and 7 . Note that if 7 


is not very close to 1, the left hand side of (76) is 0( 1), and since f3 <C 1, we must have 
n > 1. The spectrum of the parameters is then nearly continuous, as one would expect in 
the semiclassical regime. 

The value of p' is largely arbitrary, as long as it satisfies 1 < p' < /3 -1 / 2 . Once p' is 


selected, the imaginary part p" is determined by Eq. (75). And since p' 1 and (5 -C 1 it 


is easy to see from (75) that p" p'. 


With a complex momentum (73), the WDW wave function (18) has the form 


#(a,0) = 


(77) 


The corresponding probability distribution can be found in terms of the Klcin-Gordon cur¬ 
rent mm- Up to a normalization constant, it is given by 


J 




TVT*), 


(78) 
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In our minisuperspace model, the current has two components, 


J a = 


and satisfies the continuity equation 


d a J° + d^ = 0. 


With a proper normalization, the component 


J*=p'l/p(«)| 2 e- 2 " 


(79) 

(80) 


(81) 


(82) 


can be interpreted as the probability density for a at a given ’’time” 0, 


dV oc J^(a, (j))da. 


(83) 


To express the decay rate in terms of the proper time t, we find the amount A0 by which 
the field 0 increases during one period of oscillation, r ~ Using the classical equation 

of motion^ for 0, 

on « 

(84) 


2 Gp_ 
37r a 3 


and ignoring, as before, the contribution to the turning points from p, we have 


. 2Gp f 1 

A0 = —— / dt- 


AGp' r+ , i 

~- ' da - 


3tt J t ~ a(t) 3 3tt J a _ 
Expressing a from Eq. (J2| we evaluate the integral: 

1 


da(ty 


(85) 


r 

A 0 = 2 p' / 

J OL 


OL + 

da , _ 

a- V^U( 


= 2tt/3p . 


( 86 ) 


a 


We now relate the field 0 to the number of oscillations, 


N = — = 


A0 2tt/ 3p' ’ 


(87) 


so the probability in Eq. (88) becomes 


J<t> yr g-47T 0p'p"N 


i/p'WP 


( 88 ) 


Note that this is different from Eq. because of the rescalings 0 = {A-kG/Z) 1 / 2 ^ andp^ = (S/IttG) 1 / 2 ^. 
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Finally, using Eq. (75), we obtain 


J 4, oc exp (e- 2Ka N) . 


(89) 


We see that the probability for the universe to remain in the oscillating state decreases 
by a factor of e in N — e 2K ° oscillations. The characteristic lifetime of a simple harmonic 
universe is thus 


T = ‘—e 2K °, 

LU 


(90) 


with K 0 given by Eq. (71). 


VI. THE TUNNELING PROBABILITY 


In the semiclassical picture, we can think of the SHU as undergoing classical oscillations 
between the turning points a_ and a+, with some probability of tunneling through the 
barrier every time it hits the point a_. We shall now calculate this tunneling probability 
and relate it to the tunneling rate that we found in the preceding section. 

We shall focus on the case of small A, when 7 1 and the turning points are approxi¬ 

mately given by 


1 


a_ 


a+ 


(91) 

(92) 


27a 0 %nGo ’ 

27 cr 

u ] A]' 

The turning points are then widely separated, a + /a_ ~ 4y 2 1, and the form of the 

barrier between a_ and a = 0 is essentially independent of A. In this regime, we expect 
the tunneling probability to be nearly the same as for a A = 0 universe undergoing a single 
bounce at a = a_. Then the probability for SHU to remain in the oscillating phase after N 
oscillations is 

V N « (1 - Q) n « e~ QN , (93) 


where Q -C 1 is the tunneling probability for a A = 0 universe. 

In order to calculate Q, we find the semiclassical WDW wave function as we did in Sec. 
Ill, except now we only have regions I, II and III to consider. We do not need a time variable 
in this case, but the scalar field still plays a useful role of introducing a classically allowed 
region near the singularity. This allows us to impose an outgoing boundary condition at 
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a —* — oo, but we assume as before that the presence of the scalar field has little effect on 
the dynamics. 

By the same argument as in Sec. Ill, the wave function at large values of a has the form 


of (18) with f p (a ) given by Eq. (53). No boundary condition is imposed at a —> +oo, so 
we do not have any quantization condition in this case, and the momentum eigenvalue 


p can be set to be real. The wave function (53) describes an ensemble of contracting 
universes, which bounce at a = a_ and re-expand. The expanding component has a smaller 
coefficient, accounting for the fact that some universes have been lost to tunneling decay. 
The probability to avoid decay is given by the ratio of the probability fluxes for the two 


components in Eq. (53), 


1 -Q = 


.n-o 

JH<~) 


,—2K 


4 + e 


-2 K 


1 — e 


-2K 


(94) 


where left and right arrows correspond to contracting and expanding branches, respectively. 


(Note that at large a both terms in (53) are very rapidly oscillating, so any interference 


effects between the two terms become completely negligible.) Thus, we have 

Q 


-2 K 


(95) 


The tunneling exponent K can be found from Eq. (pfTj) . In the limit of large 7 it gives 


Ah 


105/3 


69 

I 67 5 


(96) 


where the second term can be dropped in the limit of 7 —y 00 . 


Substituting Q from (95) in Eq. (93), we recover Eq. (89), as expected. 


VII. CONCLUSIONS 

Our goal in this paper was to implement DeWitt’s prescription, that time evolution in 
quantum cosmology should be described in terms of semiclassical superspace variables, which 
can be used to define a “clock”. We applied this approach to the calculation of the tunneling 
decay rate of a simple harmonic universe. The role of a clock in our model was played by 
a homogeneous, massless, minimally coupled scalar field <fi. The classical evolution of (p is 
monotonic, and thus it is a good time variable. 

We found the WKB wave functions v h(a,0), which are eigenstates of the momentum p 
conjugate to (p , and matched these wave functions across the turning points, where the WKB 
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approximation breaks down. We imposed a boundary condition at a —> oo requiring that T 
vanishes in that limit and an outgoing boundary condition at a = 0. The latter condition 
means that collapse to a = 0 is irreversible, so collapsing universes do not bounce back from 
the singularity. These two boundary conditions determine the wave function completely and 
in addition provide two constraints on the parameters of the system and on the momentum 
eigenvalue p^. We showed how these constraints can be used to calculate the decay rate. 

We also considered the case of a vanishing cosmological constant A, when the universe 
experiences a single bounce off the barrier and found the tunneling probability through the 
barrier using the conserved Klcin-Gordon-type current. The resulting probability agrees 
with our calculation of the decay rate in the limit of small A. 

It should not be difficult to extend our analysis to a static universe, which has 7 = 1 
and a + = a_ = uj~ 1 . I 11 this case, the classically allowed region III reduces to a single 
point, and the method of a linear approximation for the potential U(a) around the turning 
points that we used in Sec. Ill cannot be applied. However, one can instead use a quadratic 
approximation U(a) oc (a — a *) 2 3 around the point a* = uj~ 1 . The wave function in that 
range can be expressed in terms of harmonic oscillator functions, which will then have to be 
matched to the WKB wave functions away from a*. Alternatively, it should be possible to 
fold the solution numerically. 
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